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Abstract.  Symmetric  Line  Gauss-Seidel  relaxation,  when  used  to  compute 
steady  solutions  to  the  upwind  differenced  Euler  equations  of  gas  dynamics, 
is  shown  to  be  unstable.  The  instability  occurs  for  the  long  waves.  If  SLGS 
is  used  in  a  multigrid  scheme,  stability  is  restored.  However,  the  use  of  an 
unstable  relaxation  scheme  will  not  provide  a  robust  multigrid  code.  Damped 
Symmetric  Point  Gauss-Seidel  relaxation  is  stable  and  provides  similar  multi¬ 
grid  convergence  rates  at  much  lower  cost.  However,  it  fails  if  the  flow  is  aligned 
with  the  grid  over  a  substantial  part  of  the  computational  domain.  Damped 
Alternating  Direction  Line  Jacobi  relaxation  can  overcome  this  problem. 
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1.  Introduction 

The  Euler  equations  that  describe  the  flow  of  an  inviscid  compressible  gas  can  be  integrated  in 
time  by  means  of  an  implicit  method.  This  is  desired  if  the  flow  displays  features  on  different 
scales,  or  if  the  steady  state  has  to  be  computed  efficiently.  The  implicit  formulation  gives  rise 
to  a  large  sparse  system  of  linear  equations,  which  can  be  solved  by  factorisation  methods  such 
as  the  Alternating  Direction  Implicit  method  [4]  and  Approximate  Factorisation  [1].  However, 
if  the  spatial  discretisation  is  obtained  by  upwind  differencing,  the  implicit  system  is  diagonally 
dominant  and  can  be  solved  more  efficiently  by  classical  relaxation  methods.  This  was  pointed 
out  and  explored  by  van  Leer  and  myself  [12],  and,  independently,  by  Chakravarty  [3]. 

In  [12]  we  found  that  one  particular  relaxation  method,  Symmetric  Line  Gauss-Seidel  (SLGS), 
did  not  converge  as  fast  as  expected,  and  sometimes  did  not  converge  at  all.  This  was  thought 
to  be  caused  by  the  numerical  boundary  conditions.  The  same  problem  was  encountered  in  [10]. 

Here  it  is  shown  that  the  convergence  problem  is  due  to  the  intrinsic  instability  of  the  SLGS 
scheme.  In  §2,  the  classical  von  Neumann  stability  analysis  is  carried  out  on  the  system  of 
linearised  Euler  equations  in  two  dimensions.  As  the  Fourier  modes,  used  in  the  analysis,  are  not 
the  proper  eigenfunctions  for  Gauss-Seidel  relaxation,  some  numerical  experiments  were  carried 
out.  The  results  are  presented  in  §3. 

The  instability  occurs  for  the  long  waves.  If  SLGS  is  applied  as  a  relaxation  scheme  in  a 
multigrid  code,  the  corrections  from  coarser  grids  are  sufficient  to  suppress  the  instability.  This 
is  shown  in  §4.  A  multigrid  scheme  based  on  SLGS  relaxation  has  already  been  used  in  [8]. 
However,  damped  Symmetric  Point  Gauss-Seidel  provides  similar  convergence  rates  at  a  lower 
cost,  and  damped  Alternating  Direction  Line  Jacobi  has  much  better  convergence  factors. 

The  main  results  are  summarised  in  §5. 

*  This  work  has  been  eupported  by  the  Center  for  Large  Scale  Scientific  Computing  (CLaSSiC)  Project  at 
Stanford  under  the  Office  of  Naval  Rerearch  Contract  N00014-83-K-0335. 


2.  Stability  analysis 

The  flow  of  an  ideal  inviscid  compressible  gas  can  be  described  by  the  Euler  equations.  In 
conservation  form,  these  are  given  by 


dw  df  dg 
dt  +  dx  By 


=  0. 


The  vector  of  states  w  and  the  fluxes  /  and  g  are 


(2.1a) 


(2.16) 


The  density  of  the  gas  is  denoted  by  p.  The  x-  and  y- component  of  the  velocity  are  u  and  v, 
respectively.  The  energy  E,  total  enthalpy  H,  pressure  p,  and  sound  speed  e  are  related  by 


E  =  £  +  ^(u3  +  v3),  H  =  E  +  p/p,  e2  =  yp/p. 

A  non-conservative  form  of  (2.1)  is 


(2.2) 


Here  S  =  log(p/p7)  is  the  specific  entropy. 

The  stability  analysis  will  be  carried  out  for  the  discretisation  of  the  linear  residual  operator 


L 


Adx+Bdy' 


(2.4) 


with  constant  coefficients  and  periodic  boundary  conditions.  The  operator  is  discretised  in  space 
by  upwind  differencing.  This  is  accomplished  as  follows.  The  matrix  A  is  diagonalised  by  Qi, 
according  to 


A  — Q|A1Q11, 


/  1  0  0  1\ 
[0-1  0  0  I 
1-1  001 
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For  B  we  have: 
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(2.56) 


To  accomplish  the  upwind  differencing,  the  matrix  A*  (k  =  1,2)  is  split  into  and  AJ ,  which 
contain  the  positive  and  negative  elements  of  At,  respectively.  This  implies 

A+  +  A;=At,  A+-A*=|At|.  (2.6) 
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Now  define 

A^aQ.AfQT1.  B±=Q7\±Qii: 

(2.7) 

It  follows  that 

A  =  A+  +  A~ ,  |A|  =  <J1|A1|Qr1  =  ^+-^": 

(2.8) 

- 

B  =  B+  +  B~ ,  |B|  =  (3a|Aa|gj 1  =  B+-B~. 

The  discrete  linear  residual  operator  becomes 


Lk  =  ^[A+(l  -  771)  +  A~(T,  -  1)]  +  i[B+(l  -  T-1)  +  B~(Ty  -  1)].  (2.9) 

Here  the  shift  operators  are  defined  by  Txwkukj  =  u>i,+i,*j*  Ttwkl  ki  =  For  simplicity, 

the  grid  is  assumed  to  be  uniform  (hx  =  h9  =  h). 

For  the  stability  analysis  we  consider  the  usual  Fourier  modes  of  the  form 

exp  [-i  (M«  +  *a0y)] ,  (2.10a) 

where  hi  and  k2  are  spatial  indices  on  a  Ni  x  N2  grid.  The  frequencies  for  a  grid  of  the  same  size 
are 

9,  =  2*^-,  9,  =  2x-£-t  hmO,...,Ni-l,  /2  =  0 . ATj-1.  (2.106) 

The  Fourier  transforms  of  the  shift-operators  Tx  and  Tt  are 

ta  =  exp(»0r),  Ty  =  exp (i0f),  0  <  9X  <  2x,  0  <  0t  <  2 zr.  (2.11) 

The  relaxation  operator  for  Symmetric  Line  Gauss-Seidel,  with  the  line  in  the  y-direction,  is 
the  product  of  a  forward  sweep 

61  =  /  -  Mkl‘Lh,  Mi  =  Lk  -  ^A~fx,  (2.12  a) 

and  a  backward  sweep 

67  -  /  -  M2lLk,  Ma  =  lk  +  (2.126) 

resulting  in 

dSLOS  zs  d2di.  (2.12c) 

Here  /  denotes  a  4  x  4  identity  matrix.  If  Mi  or  M2  is  singular,  the  pseudo-inverse  should  be 
used. 

The  stability  of  this  scheme  is  investigated  by  computing  the  amplification  factor  of  the 
residual 

i^smax^.f,),  xr(^.^)sp(L',G5£G5(L'k)t).  (2.13) 

Here  p(  )  denotes  the  spectral  radius.  The  operator  Lk  can  be  singular  [9:  Lemma  3.1],  and  the 
waves  for  which  it  is  singular  are  obviously  not  damped  or  amplified.  To  exclude  these  waves, 
the  operator  Lk  and  its  pseudo-inverse  are  included  in  the  definition  of  kt  ■  The  value  of  7cr  will 
depend  on  the  velocities  u/c  and  v/c,  and  on  the  grid-size  Nx  x  N2. 

SLGS  is  an  exact  solver  for  the  fourth  equation  of  (2.3),  which  describes  the  convection  of 
entropy  along  streamlines.  It  is  also  exact  if  |u/e|  >  1. 
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Figure  1  shows  Kr($a,0t)  for  u/c  —  0.8  and  v/e  =  0.  Here  N\  =  Ni  =  32.  The  instability  is 
clearly  visible.  It  occurs  for  the  long  waves,  at  \6X\  =  |0V|  =  2x/32.  For  the  given  u/c  and  v/e, 
the  instability  does  not  show  up  on  a  16  x  16  grid,  whereas  it  becomes  worse  on  finer  grids. 

Similar  instabilities  occur  for  most  values  of  |u/c|  <  1  and  |v/c|  <  1,  in  some  cases  on  grids 
coarser  than  in  the  example  of  Fig.  1. 

3.  Numerical  experiments 

It  is  well  known  that  Fourier  modes  are  not  the  proper  eigenfunctions  for  Gauss-Seidel  relaxation. 
Therefore,  the  validity  of  the  above  analysis  may  be  questioned,  especially  for  the  longer  waves. 
In  this  section  the  instability  will  be  investigated  by  numerical  experiments  on  the  system  of 
nonlinear  Euler  equations  (2.1). 

For  the  upwind  differencing,  van  Leer’s  flux-vector  splitting  (FVS)  [11]  is  used  as  an  ap¬ 
proximate  Riemann-solver.  This  scheme  gives  rise  to  matrices  A*  =  df^/dw  and  B±  =  dg±/dw 
which  are  different  from  those  in  (2.7).  Therefore,  the  stability  properties  of  this  scheme  will  be 
different  from  those  predicted  in  the  previous  section,  but  not  by  too  much.  Figure  2  shows  the 
amplification  factor  for  FVS,  using  the  same  parameters  as  in  Fig.  1. 

As  a  test  problem,  flow  through  a  straight  channel  is  considered.  The  grid  is  square  and 
uniform.  There  are  hard  walls  on  the  lower  and  upper  side.  Boundary  conditions  at  the  walls 
are  implemented  by  ghost  cells  that  contain  reflected  states.  Characteristic  boundary  conditions 
are  used  at  the  inlet  and  outlet.  In  principle,  overspecification  can  be  used  because  the  Riemann- 
solver  takes  care  of  the  appropriate  switching  between  incoming  and  outgoing  characteristics. 
However,  because  FVS  is  not  a  very  good  approximate  Riemann-solver,  the  use  of  characteristic 
boundary  conditions  is  recommended. 

The  free-stream  values  are  chosen  to  be 

Poo  —  1>  Uoo  =  0*8,  Woo  —  0,  Coo  ~  1>  (8-1) 

whereas  the  gas-constant  7  =  1.4.  A a  initial  conditions,  we  take  the  free-stream  values  and  add 
random  noise  with  an  amplitude  of  10~s.  The  steady  state  is  given  by  the  free-stream  values. 

Table  1  shows  the  predicted  and  observed  amplification  factors.  There  is  a  clear  qualitative 
agreement.  Inspection  of  the  error,  which  is  defined  at  the  difference  between  the  non-converged 
solution  and  the  numerical  steady  state,  confirms  the  long-wave  instability. 


Table  1.  Amplification  factors  for  flux- vector  splitting  on 
grids  of  various  sizes,  with  Uoo/coo  =  0.8  and  Voo/coo  = 
0.0.  The  result  of  the  Fourier  stability  analysis  is  denoted 
by  Kr.  T he  other  values  are  determined  from  numerical 
experiments  on  the  full  system  of  nonlinear  Euler  equations. 


4.  Multigrid 

The  instability  of  Symmetric  Line  Gauss-Seidel  occurs  for  the  long  waves.  It  is  therefore  expected 
that  the  combination  of  SLGS  and  the  multigrid  technique  will  provide  a  stable  scheme.  The 
analysis  of  multigrid  convergence  for  the  linearised  Euler  equations  with  constant  coefficients  is 
presented  in  detail  in  [9].  The  multigrid  convergence  factor  XP  of  a  given  relaxation  scheme  is 
estimated  by  considering  two  grids,  a  fine  and  a  coarse.  The  number  of  cells  on  the  coarse  grid  is 
1/4  of  that  on  the  fine.  For  the  restriction  to  the  coarser  grid,  simple  averaging  is  used.  Zero-order 
interpolation  is  applied  for  the  prolongation  back  to  the  fine  grid  (cf.[7]).  In  the  analysis,  it  is 
assumed  that  the  coarse-grid  equations  are  solved  exactly.  The  combined  result  of  the  coarse-grid 
correction  and  the  relaxation  scheme  is  described  by  Xr.  This  two-level  multigrid  convergence 
factor  gives  a  reasonable  estimate  of  the  convergence  speed  when  many  coarser  grids  are  used. 

The  result  of  the  so-called  two-level  analysis  is  shown  in  Fig.  3.  The  instability  has  disap¬ 
peared.  The  overall  convergence  rate  is  good,  except  near  the  singularities  of  the  (linear)  residual. 
The  bad  convergence  around  v  =  0  occurs  when  the  flow  is  aligned  with  the  grid  over  a  substantial 
part  of  the  computational  domain.  To  overcome  the  problem  of  strong  alignment  [2],  one  might 
consider  a  relaxation  scheme  that  consists  of  a  SLGS  step  with  the  line  in  the  z-direction,  followed 
by  SLGS  with  the  line  in  the  p-direction.  This  will  be  called  Alternating  Direction  Symmetric 
Line  Gauss-Seidel.  ADSLGS  turns  out  to  be  stable  in  a  larger  region  of  the  (u/c,  v/c)-plane. 
However,  the  instability  persists  for  some  parameters  and  is  so  severe  that  it  does  not  disappear 
in  a  multigrid  scheme. 

A  scheme  that  suffers  from  strong  alignment  in  a  similar  way  as  SLGS,  is  damped  Symmetric 
Point  Gauss-Seidel  (SGS)  [9].  It  is  stable  as  a  single-grid  scheme  and  much  cheaper  than  its 
undamped  line  variant,  and  is  therefore  to  be  preferred.  If  one  wants  a  uniformly  good  convergence 
rate,  a  multigrid  scheme  that  employs  damped  Alternating  Direction  Line  Jacobi  (ADLJ)  can 
be  used.  The  linear  two-level  analysis  predicts  a  multigrid  convergence  rate  Xr(u/c,  v/c)  <  0.526 
for  this  scheme.  The  relaxation  operator  is  described  by  the  product  of  line  relaxation  in  the 
p-direction,  namely, 

Gx  =  f-  H;lLh,  Hi  —  Lh  +  +  1)  -  A~(Tt  +  1)],  (4.1a) 

and  line  relaxation  in  the  z-direction: 

Oi  =  I-HjlLh,  H7  =  Ll,  +  ^[B+(f-1  +  l)-B-(f,  +  l)].  (4.16) 

Note  that  the  relaxation  operators  Hi  and  H7  are  obtained  by  selecting  the  main-diagonal  (in 
terms  of  blocks)  and  two  off-diagonals  (in  one  direction)  from  the  residual  operator.  The  damping 
is  obtained  by  subtracting  the  two  other  off-diagonals  in  the  direction  perpendicular  to  the  line 
from  the  main-diagonal.  Figure  4  shows  the  two-level  multigrid  convergence  rate  for  damped 
ADLJ  relaxation. 

5.  Conclusions 

Symmetric  Line  Gauss-Seidel  relaxation  is  unstable  and  should  therefore  not  be  used  for  the 
implicit  Euler  equations.  The  instability  occurs  for  the  long  waves.  It  disappears  in  a  multigrid 
scheme.  However,  in  that  case  damped  Symmetric  Point  Gauss-Seidel  is  a  better  choice,  because 
it  is  stable  as  a  single-grid  scheme  and  can  produce  similar  multigrid  convergence  factors  at  a 
much  lower  cost  [9].  This  scheme,  however,  does  not  provide  good  convergence  in  cases  of  strong 
alignment,  the  flow  being  aligned  with  the  grid. 

A  uniformly  good  convergence  rate  can  be  obtained  with  damped  Alternating  Direction  Line 
Jacobi,  which  has  about  the  same  cost  as  SLGS.  The  two-level  analysis  indicates  a  convergence 
factor  that  is  at  worst  0.562. 


The  use  of  SLGS  has  been  recommended  for  the  Navier-Stokes  equations  in  [5].  Here  the 
”  Euler  part”  of  the  equations  is  discretised  by  upwind  differencing.  Because  SLGS  is  unstable  for 
the  long  waves,  it  is  expected  that  the  additional  viscous  terms  will  not  stabilise  the  scheme.  This 
casts  serious  doubts  on  the  general  validity  of  the  grid-independent  convergence  rates  claimed  in 
[61- 
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Fig.  3.  Multigrid  convergence  factor  for  SLGS  relaxation  aa  a  function  of  u  and  v  (c  —  l).  Each 
point  is  computed  for  a  64  x  64  grid. 


Fig.  4.  Multigrid  convergence  factor  for  ADLJ  relaxation  as  a  function  ofu  and  v  (c  =  1 ).  Each 
point  is  computed  for  64  x  64  grid. 


tv 

;y 


v< 

M- 

$ 

V* 

k 


6 


Kl. 

Kt 


I, 


•1 

i 


L 


